Probability distribution of returns in the Heston model with stochastic volatility 
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We study the Heston model, where the stock price dynamics is governed by a geometrical (mul- 
tiplicative) Brownian motion with stochastic variance. We solve the corresponding Fokker-Planck 
equation exactly and, after integrating out the variance, find an analytic formula for the time- 
dependent probability distribution of stock price changes (returns). The formula is in excellent 
agreement with the Dow- Jones index for the time lags from 1 to 250 trading days. For large returns, 
the distribution is exponential in log-returns with a time-dependent exponent, whereas for small 
returns it is Gaussian. For time lags longer than the relaxation time of variance, the probability 
distribution can be expressed in a scaling form using a Bessel function. The Dow- Jones data for 
1982-2001 follow the scaling function for seven orders of magnitude. 
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I. INTRODUCTION 

Stochastic dynamics of stock prices is commonly de- 
scribed by a geometric (multiplicative) Brownian motion, 
which gives a log-normal probability distribution func- 
tion (PDF) for stock price changes (returns) How- 
ever, numerous observations show that the tails of the 
PDF decay slower than the log-normal distribution pre- 
dicts (the so-called "fat-tails" effect) |§ § §. Particu- 
larly, much attention was devoted to the power-law tails 
^ M. The geometric Brownian motion model has two 
parameters: the drift fx, which characterizes the aver- 
age growth rate, and the volatility a, which characterizes 
the noisiness of the process. There is empirical evidence 
and a set of stylized facts indicating that volatility, in- 
stead of being a constant parameter, is driven by a mean- 
reverting stochastic process |Q |j| . Various mathematical 
models with stochastic volatility have been discussed in 
literature § @ || g| || 0, Q. 

In this paper, we study a particular stochastic volatil- 
ity model, called the Heston model Q , where the square 
of the stock-price volatility, called the variance v, follows 
a random process known in financial literature as the 
Cox-Ingersoll-Ross process and in mathematical statis- 
tics as the Feller proces s [p] , fL6| . Using the Fourier and 
Laplace transforms jli], [l6f , we solve the Fokker-Planck 
equation for this model exactly and find the joint PDF 
of returns and variance as a function of time, conditional 
on the initial value of variance. While returns are read- 
ily known from a financial time-series data, variance is 
not given directly, so it acts as a hidden stochastic vari- 
able. Thus, we integrate the joint PDF over variance and 
obtain the marginal probability distribution function of 
returns unconditional on variance. The latter PDF can 
be directly compared with financial data. We find an 
excellent agreement between our results and the Dow- 
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Jones data for the 20-years period of 1982-2001. Us- 
ing only four fitting parameters, our equations very well 
reproduce the PDF of returns for time lags between 1 
and 250 trading days. In contrast, in ARCH, GARCH, 
EGARCH, TARCH, and similar models, the number of 
fitting parameters can easily go to a few tens [JlTj . 

Our result for the PDF of returns has the form of a 
one-dimensional Fourier integral, which is easily calcu- 
lated numerically or, in certain asymptotical limits, an- 
alytically. For large returns, we find that the PDF is 
exponential in log-returns, which implies a power-law dis- 
tribution for returns, and we calculate the time depen- 
dence of the corresponding exponents. In the limit of 
long times, the PDF exhibits scaling, i.e. it becomes a 
function of a single combination of return and time, with 
the scaling function expressed in terms of a Bessel func- 
tion. The Dow- Jones data follow the predicted scaling 
function for seven orders of magnitude. 

The original paper |pd| solved the problem of option 
pricing for the Heston model. Numerous subsequent 
studies [|l3| [l4|, [l5[ [l^] compared option pricing derived 
from this model and its extensions with the empirical 
data on option pricing. They found that the Heston 
model describes the empirical option prices much better 
than the Black-Scholes theory, and modifications of the 
Heston model, such as adding discontinuous jumps, fur- 
ther improve the agreement. However, these papers did 
not address the fundamental question whether the stock 
market actually follows the Heston stochastic process or 
not. Obviously, if the answer is negative, then using the 
Heston model for option pricing would not make much 
sense. The stock-market time series was studied in Ref. 
p5| jointly with option prices, but the focus was just on 
extracting the effective parameters of the Heston model. 
In contrast, we present a comprehensive comparison of 
the stock market returns distribution with the predic- 
tions of the Heston model. Using a single set of four 
parameters, we fit the whole family of PDF curves for 
a wide variety of time lags. In order to keep the model 
as simple as possible with the minimal number of fitting 
parameters, we use the original Heston model and do not 
include later modifications proposed in literature, such as 
jumps, multiple relaxation time, etc. [^3[ |lj, |lq| . Inter- 



2 



estingly, the parameters of the model that we find from 
our fits of the stock market data are of the same order of 
magnitude as the parameters extracted from the fits of 
option prices in Refs. [l4[ 



II. THE MODEL 

We consider a stock, whose price St, as a function of 
time t, obeys the stochastic differential equation of a geo- 
metric (multiplicative) Brownian motion in the Ito form 
11,01: 



dS t = pS t dt + a t S t dW t 



(i) 



(1) 



Here the subscript t indicates time dependence, p is the 
drift parameter, W t is a standard random Wiener pro- 
cess 1 , and at is the time-dependent volatility. 

Since any solution of (^) depends only on of, it is con- 
venient to introduce the new variable Vt = of, which is 
called the variance. We assume that Vt obeys the follow- 
ing mean-reverting stochastic differential equation: 



dvt = -j(v t -0)dt + K*Jv~tdW K t 



(2) 



(2) 



Here 6 is the long-time mean of v, 7 is the rate of relax- 
(2) 

ation to this mean, W t is a standard Wiener process, 
and k is a parameter that we call the variance noise. Eq. 
(||) is known in financial literature as the Cox-Ingersoll- 
Ross (CIR) process and in mathematical statistics as the 
Feller process H, Rq. Alternative equations for v±, with 

the last term in (0) replaced by k dW^ or nv t dW^ , 
have been also discussed in literature However, in 
our paper, we study only the case given by Eq. (2|). 

We take the Wiener process appearing in (p) to be 
correlated with the Wiener process in (0): 



dW} 



pdW t {1) + yJ\-p 2 dZ u 



(3) 



where Z t is a Wiener process independent of Wt , and 
p e [— 1, 1] is the correlation coefficient. A negative cor- 



relation (p < 0) between and W t w is known as the 
leverage effect ||, p. 41]. 

It is convenient to change the variable in (Q) from price 
St to log-return r t = ln(S t /So). Using Ito's formula [^9| , 
we obtain the equation satisfied by r t : 



r(2) 



dr. 



(ji-^dt + yfvldwP. 



(4) 



The parameter p can be eliminated from (^) by changing 
the variable to x t = r t — pt, which measures log-returns 
relative to the growth rate p: 



d Xt = -^dt + ^r t dw t {1) . 



(5) 



1 The infinitesimal increments of the Wiener process dWt are 
normally-distributed (Gaussian) random variables with zero 
mean and the variance equal to dt. 




FIG. 1: The stationary probability distribution Ti*(v) of vari- 
ance v, given by Eq. (g) and shown for a — 1.3 from Table [l| 
The vertical line indicates the average value of v. Inset: The 
corresponding stationary probability distribution ill <T '(ii) of 
volatility a given by Eq. (|l0|). 



Where it does not cause confusion with r t , we use the 
term "log-return" also for the variable Xt- 

Equations (||) and (g) define a two-dimensional 
stochastic process for the variables x t and v t jll], pT[ . 
This process is characterized by the transition probabil- 
ity P t (x,v \vi) to have log-return x and variance v at 
time t given the initial log-return x = and variance Uj 
at t = 0. Time evolution of Pt {x, v \ m) is governed by the 
Fokker-Planck (or forward Kolmogorov) equation [fill] 



^=^i(..-w + i|V>> 
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+PK 



d 2 



dx dv 



(vP) 



1 d 2 

2dx 2 



(vP) 



2 dv 2 



(vP). 



The initial condition for (||) is a product of two delta 
functions 



Pt=o(x,v\ Vi) = 5(x) S(v - Vj). 



(7) 



The probability distribution of the variance itself, 
n t (w) = J dx Pt(x, v ), satisfies the equation 



0_ 

dt 



u t (v) = |- [ 7 (« - e)u t (v)} + i vU ^ . ( 8 ) 



which is obtained from (Q) by integration over x. Feller 
p6[ has shown that this equation is well-defined on the 
interval v E [0, +00) as long as 9 > 0. Eq. (||) has the 
stationary solution 



n.(«) 



T(a) 



,-av/6 



«=% (9) 



which is the Gamma distribution. The parameter a in 
(^) is the ratio of the average variance 9 to the character- 
istic fluctuation of variance k 2 /2 r y during the relaxation 
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time I/7. When a — > 00, II* (u) — > 5(u — 8). The corre- 
sponding stationary PDF of volatility a is 



y ' T(a) 9 a 



(10) 



Functions (|9j) and (|10|) are integrable as long as a > 0. 

The distributions II* (u) and Ili CT ' > (cr) are shown in Fig. |l| 
for the value a = 1.3 deduced from the fit of the Dow- 



Jones time series and given in Table m in Sec. VIII. 



III. SOLUTION OF THE FOKKER-PLANCK 
EQUATION 

Since x appears in (|^) only in the derivative operator 
d/dx, it is convenient to take the Fourier transform 

P t (x,v\v i )= ^e^ s P t)P >k). (11) 

J-00 27r 

Inserting ([ll]) into (^), we find 

y pl — ip x . d k 2 d 2 



(12) 



Eq. (|12|) is simpler than (|6j) , because the number of vari- 
ables has been reduced to two, v and i, whereas p x only 
plays the role of a parameter. 

Since Eq. ( |l2] ) is linear in v and quadratic in d/dv, it 
can be simplified by taking the Laplace transform over v 



Pt,p x (Pv I Vi) 



dve- p " v Pt, Px (v\vi). 



(13) 



The partial differential equation satisfied by Pt,p* {p v \ Vi) 
is of the first order 



d_ 

Of 



Fp v + 



Px - l Px 



d 



P = -j9 Pv P, 
(14) 



2 J dp 

where we introduced the notation 

T = j + ipKp x . (15) 
Eq. ( |l4| ) has to be solved with the initial condition 

-P*=o, Px (Pv I Vi) = exp(-p t ,w i ). (16) 

The solution of @ is given by the method of charac- 
teristics ppf : 

P^pAPv I Vi) = exp (-p v (0)v t - 76 J dTp v (T)^j , (17) 

where the function p v (r) is the solution of the character- 
istic (ordinary) differential equation 

^=r^(r) + ^p1(r)-^ (18) 



with the boundary condition p v (t) — p v specified at r = 
t. The differential equation (|18|) is of the Riccati type 
with constant coefficients [Elf , and its solution is 



Pv(t) 



2Q 



1 



k 2 Ce°(*- T ) - 1 
where we introduced the frequency 



T-fl 



n = ^v 2 + k 2 {p 2 x 

and the coefficient 

c = i 2n 



lPx) 



k 2 Pv + (r - n) ' 

Substituting @ into @, we find 

Pt.pSPv I «i) 

= exp^ -p„(0)^+ r 1 ' Wn- 



(19) 



(20) 



(21) 



(22) 



-m 



c-i 



IV. AVERAGING OVER VARIANCE 

Normally we are interested only in log-returns x and 
do not care about variance v. Moreover, whereas log- 
returns are directly known from financial data, variance 
is a hidden stochastic variable that has to be estimated. 
Inevitably, such an estimation is done with some degree 
of uncertainty, which precludes a clear-cut direct com- 
parison between Pt{x, v \ Vi) and financial data. Thus we 
introduce the reduced probability distribution 

+00 

Pt(x\vi)= fdvP t (x,v\vi) = I ^e^iv(OK), 


(23) 

where the hidden variable v is integrated out, so p v = 0. 
Substituting £ from ( |2l| ) with p v — into (j2^) , we find 

P t (x\Vi)= [ — e iVxX ~ Vi r+o *otMm/2) 



x e 



2tt 



2^ ln(cosh qi + g sinh ^) + ^ 



(24) 



To check the validity of ( p24j ) , let us consider the limit- 
ing case k = 0. In this case, the stochastic term in (^) is 
absent, so the time evolution of variance is deterministic: 

v t = 6+(vi- 9)6-1*. (25) 

Then process (||) gives a Gaussian distribution for x, 



P 



(k=0) 



(i£ | Vi) = 



V27rtf 



= ex P 



(x + Vtt/2) 2 
2v t t 



(26) 



with the time-averaged variance v t = \ J * dr v T . On the 
other hand, by taking the limit k — > and integrating 
over p x in (|24|), we reproduce the same expression (|26|). 
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Dow Jones data, 1982-2001 and 1990-2001 
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FIG. 2: Probability distribution Pt(x) of log-return x for 
different time lags t. Points: The 1982-2001 Dow- Jones data 
for t = 1, 5, 20, 40, and 250 trading days. Solid lines: Fit of 
the data with Eqs. ( pfi| ) and (p9|). For clarity, the data points 
and the curves for successive t are shifted up by the factor 
of 10 each. Inset: The 1990-2001 Dow- Jones data points 
compared with the same theoretical curves. 



Eq. ( p4| ) cannot be directly compared with financial 
time-series data, because it depends on the unknown ini- 
tial variance Vi. In order to resolve this problem, we 
assume that Vi has the stationary probability distribu- 
tion II* (•»,), which is given by (g). Thus we introduce 
the probability distribution function Pt(x) by averaging 
(B4T) over Vi with the weight II* (uj): 



Pt{x) 



dvi U*(vi) Pt(x | Vi). 



(27) 



The integral over v\ is similar to the one of the Gamma 
function and can be taken explicitly. The final result is 
the Fourier integral 



Pt(x) 



1 

2^ 



dp x e 



ipxX+Ft(p x ) 



with 



F t ( Px ) = I 1 rt 



(28) 



(29) 



2 7 6» 



In 



nt n 2 - r 2 + 2 7 r , ot 

cosh — H — — smh — 

2 2yfl 2 



The variable p x enters (p9| ) via the variables T from 
(H) and Q from ©. It is easy to check that Pt(x) is 
real, because Re_F is an even function of p x and ImF is 
an odd one. One can also check that Ft(p x = 0) = 0, 
which implies that Pt{x) is correctly normalized at all 
times: J dxPtix) — 1. The simplified version of Eq. j29| ) 
for the case p — isgiven in Appendix |A|. 

Eqs. ( pjj ) and (^9|) for the probability distribution 
Pt(x) of log-return x at time t are the central analyti- 



cal result of the paper. The integral in ( |28| ) can be cal- 
culated numerically or, in certain regimes discussed in 
Sees. |y|, [v|, and VII , analytically. In Fig. ||, the calcu- 
lated function Pt{x), shown by solid lines, is compared 
with the Dow-Jones data, shown by dots. (T echnic al de- 
tails of the data analysis are discussed in Sec. VIII .) Fig. 
H demonstrates that, with a fixed set of the parameters 
7, 8, k, [i, and p, Eqs. ( p8|) and (p9] ) very well reproduce 
the distribution of log-returns x of the Dow- Jones index 
for all times t. In the log-linear scale of Fig. ||, the tails 
of hiP t (x) vs. x are straight lines, which means that that 
tails of Pt{x) are exponential in x. For short times t, 
the distribution is narrow, and the slopes of the tails are 
nearly vertical. As the time t progresses, the distribution 
broadens and flattens. 



V. ASYMPTOTIC BEHAVIOR FOR LONG 
TIME t 

Eq. (j^) implies that variance reverts to the equilibrium 
value 9 within the characteristic relaxation time 1/7. In 
this section, we consider the asymptotic limit where time 
t is much longer than the relaxation time: jt ^> 2. Ac- 
cording to ( J15| ) and (p0| ), this condition also implies that 
fit > 2. Then Eq. d29T reduces to 



F t ( Px )^^(T-n). 



(30) 



Let us change of the variable of integration in ( |28|) to 



where 
Po : 



Px 



2p 7 



LOq 



■■Px +IP0, 



(31) 



7 2 + K 2 (1 _ p 2 )p 2, (32) 



Substituting (|l|) into @, @, and @, we transform 
(p8|) to the following form 



Pt{x) 

where 
A = 

and 



-p x+At roo 

; / dp x cos(Ap x )e- B V^+Pi 

1 - p A Jo 



ujoe 

TTKy/l — p 



«V1 



y6t 

, x + p 

n2 \ K 



B 



A 



7f? 27 — pK 



(33) 



(34) 



(35) 



2k 2 1-p 2 ' 

According to formula 3.914 from [f22| , the integral in (|3_ 
is equal to BATi(V^4 2 + B 2 )/^JA 2 + B 2 , where ifi is the 
first-order modified Bessel function. 

Thus, Eq. ( p8| ) in the limit 7 t 3> 2 can be represented 
in the scaling form 

P t (x) = N t e- poX P*(z), P*(z) = K x {z)jz, (36) 
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Dow-Jones data, 1982-2001 and 1990-2001 
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FIG. 3: Renormalized probability density P t {x)e p ° x /N t plot- 
ted as a function of the scaling argument z given by (37). 
Solid lines: The scaling function P*(z) = Ki(z)/z from (36), 
where K\ is the first-order modified Bessel function. Upper 
and lower sets of points: The 1982-2001 and 1990-2001 Dow- 
Jones data for different time lags t. For clarity, the lower data 
set and the curve are shifted by the factor of 10~ 2 . 
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FIG. 4: The fraction Gt of the total probability contained 
in the Gaussian part of Pt(x) vs. time lag t. Inset: Time 
dependence of the probability density at maximum Pt(x m ) 
(points), compared with the Gaussian f -1 / 2 behavior (solid 
line). 
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where the argument z = ^/ A 2 + B 2 is 



{x + fTfOt/Kf 
I-P 2 



and the time-dependent normalization factor Nt is 

L0 2 j9t 



! 7i" 



(38) 



Eq. (|3q ) demonstrates that, up to the factors Nt and 
e~ p ° x , the dependence of Pt(x) on the two arguments x 
and t is given by the function P*(z) of the single scaling 
argument z in (|37|). Thus, when plotted as a function 
of z, the data for different x and t should collapse on 
the single universal curve P*(z). This is beautifully illus- 
trated by Fig. |j| where the Dow- Jones data for different 
time lags t follows the curve P*(z) for seven orders of 
magnitude. 

In the limit z 1, we can use the asymptotic expres- 
sion [|| Ki(z) « e- z v / n/2z in @ and take the loga- 
rithm of P. Keeping only the leading term proportional 
to z and omitting the subleading term proportional to 
lnz, we find that lnP t (a;) has the hyperbolic distribution 
§ P- 14] 



In 



N 



~p a x 



for z>l. 



(39) 



Let us examine Eq. ( p9[ ) for large and small |x|. 

In the first case |xf 3> jOt/n, Eq. (37) gives z 
wo\x\/ny/l — p 2 , so Eq. (^Tj|) becomes 



In 



Nt 



-p x 



UJ 



- P 2 



\x\. 



(40) 



Thus, the PDF Pt(x) has the exponential tails (]4(]) for 
large log-returns \x\. Notice that, in the considered limit 
jt 2, the slopes d\nP/dx of the exponential tails (40) 
do not depend on time t. Because of po, the slopes (4^) 
for positive and negative x are not equal, thus the dis- 
tribution Pt(x) is not symmetric with respect to posi- 
tive and negative price changes. According to (|32|) , this 
asymmetry is enhanced by a negative correlation p < 
between stock price and variance. 

In the second case \x + p~y9t/ 'k\ -C j6t/n, by Taylor- 
expanding z in ( |37"| ) near its minimum in x and substi- 
tuting the result into (|39|), we get 



In 



Pt(x) 



-p x 



uj (x + p^9t/ny 
2(1- p 2 )j9t 



(41) 



where N{. 



N t ex-p(~uJo'y9t/K 2 ). Thus, for small log- 
returns \x\, the PDF Pt(x) is Gaussian with the width 
increasing linearly in time. The maximum of Pt(x) in 
([ll]) is achieved at 



Xm{t) = "2^ 1 1 



P(U} - 7) 



(42) 



Eq. (^2) gives the most probable log- return x m (t) at time 
t, and the coefficient in front of t constitutes a correction 
to the average growth rate p, so that the actual growth 
rate is p = p + dx m /dt. 

As Fig. g illustrates, \nP t (x) is indeed linear in x for 
large |x| and quadratic for small \x\, in agreement with 
(p0|) and (pi]). As time progresses, the distribution, which 
has the scaling form ( prj] ) and (|37|), broadens. Thus, 
the fraction G t of the total probability contained in the 
parabolic (Gaussian) portion of the curve increases, as 
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FIG. 5: Complex plane of p x . Dots: The singularities of 
Ft{p x )- Circled crosses: The accumulation points ±iq^r of 
the singularities in the limit yt 2> 2. Cross: The saddle point 
p s , which is located in the upper half-plane for x > 0. Dashed 
line: The contour of integration displaced from the real axis 
in order to pass through the saddle point p 3 . 

illustrated in Fig. ||. (The procedure of calculating G t is 
explained in Appendix [b].) Fig. |] shows that, at suffi- 
ciently long times, the total probability contained in the 
non-Gaussian tails becomes negligible, which is known 
in literature j^] . The inset in Fig. || illustrates that the 
time dependence of the probability density at maximum, 
Pt(xm), is close to i -1 / 2 , which is characteristic for a 
Gaussian evolution. 



VI. ASYMPTOTIC BEHAVIOR FOR LARGE 
LOG-RETURN x 

In the complex plane of p x , function F(p x ) becomes 
singular at the points p x where the argument of the log- 
arithm in ( p9| ) vanishes. These points are located on the 
imaginary axis of p x and are shown by dots in Fig. |^. 
The singularity closest to the real axis is located on the 
positive (negative) imaginary axis at the point pf (p^)- 
Because the argument of the logarithm in ( |29| ) vanishes at 
these two points, we can approximate F(p x ) by the dom- 
inant, singular term: F(p x ) « — (2^9/k 2 ) ln(p x — pf). 

For large |x|, the integrand of fl28[ ) oscillates very fast 
as a function of p x . Thus, we can evaluate the integral 
using the method of stationary phase |H | by shifting the 
contour of integration so that is passes through a saddle 
point of the argument ip x x + F(p x ) of the exponent in 
( p8| ) . The saddle point position p s , shown in Fig. |s| by 
the cross, is determined by the equation 



dF(p x ) 



dp x 



2 7 6> 
ft- 2 



x 



— L r, x > 0, 
Ps-vt (43) 



For a large \x\, such that \xpf\ 2j6/k , the saddle 



l. Ps~P 1 

±1 ^ /.,--> 



x < 0. 



FIG. 6: Solid line: The slope qf = -dlnP/dx of the 
exponential tail for x > as a function of time. Points: 
The asymptotic approximation ( fjrj| ) for the slope in the limit 
yt <C 2. Dashed line: The saturation value qt for jt ^> 2, 
Eq. ©. 



point p s is very close to the singularity point: p s p^ 
for x > and p s ~ p^ for x < 0. Then the asymptotic 
expression for the probability distribution is 



Ptix) 



e- xq t , x > Q, 
e xq t , x < 0, 



(44) 



where qf — ^fipf(t) are real and positive. Eq. ( fill ) shows 
that, for all times t, the tails of the probability distribu- 
tion Pt(x) for large |x| are exponential. The slopes of 
the exponential tails, — =F d In P/dx, are determined 
by the positions pf of the singularities closest to the real 
axis. 

These positions pf{t) and, thus, the slopes qf 1 depend 
on time t. For times much shorter than the relaxation 
time (yt <C 2), the singularities lie far away from the real 
axis. As time increases, the singularities move along the 
imaginary axis toward the real axis. Finally, for times 
much longer than the relaxation time (yt ^> 2), the sin- 
gularities approach limiting points: pf — > ±iq^, which 
are shown in Fig. || by circled crosses. Thus, as illus- 
trated in Fig. [| the slopes qf monotonously decrease in 
time and saturate at long times: 



1/ 



±p for 7i>2. 



(45) 



The slopes (ft5|) are in agreement with Eq. (M) valid for 



jt 3> 2. The time dependence q t 
also found analytically: 



at short times can be 



k y t 



(46) 



The dotted curve in Fig. || shows that Eq. ( [46"| ) works 
very well for short times t, where the slope diverges at 
t -> 0. 
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VII. ASYMPTOTIC BEHAVIOR FOR SHORT 
TIME t 

For a short time t, we expand the equations of Sec. 
[II to the first order in t and set p v = 0. The last term 
in Eq. ( p2| ) cancels the penultimate term, and Eq. JI§| ) 
gives Pv(0) — t(p x — ipx)/2. Substituting this formula 
into (E2f) and taking the integral (E3|) over p x , we find 



Pt(x\vi) = 



1 



y/2-KVit 

Eq. (|47| ) shows that, for a short t, the probability distri- 
bution of x evolves in a Gaussian manner with the initial 
variance Vi, because variance has no time to change. 
Substituting @ and (§) into (J27j), we find 

a „-x/2 poo 

P t (x)= % . g= = / dfrff-V^'-*/*' (48) 

where 5, = Vj/0, A = a + 6t/8, B = x 2 /29t, and C = 
a - 1/2. According to formula 3.471.9 from j22), the 
integral in © is 2(P/A) C / 2 X C (2\/AB) for RcA > 
and Rc£? > 0, where Kc is the modified Bessel function 
of the order C. Taking into account that A ss a (because 
t <C I67/K 2 for short t), we obtain the final expression 



Pt{x) 



)1 — a , 




r(a) V 

where we introduced the scaling variable 

2ax 2 _ 2^7 |x| 
6t ~ ~K~ ~ji 



a - 1/2 K a _ 1/2 (y), (49) 



v 



(50) 



In the limit y 3> 1, using the formula K u (y) 
e-yy/n/2y in @, we find 



Pt{x) 



nl/2-a 

r(«) 



0i y 



(51) 



Eqs. © and © show that the tails of the distribu- 
tion are exponential in x, and the slopes dhiP/dx are in 
agreement with Eq. (fl6|). 

In the opposite limit 1/ < 1, the small argument ex- 
pansion of the Bessel function can be found from the 
following equations pSj: 



and 



7T I- V {y) - h(y) 
2 sin(^7r) 



sin(7ri/) 



= i»r(i - u), 

(52) 



1 Ay) 



(2/ 2 /4) fc 



^ k\T(v + k+ 1)' 



(53) 



Substituting (|53| ) into (p2|), we find in the case 1/2 < a < 
3/2 



r(a-l/2) /2/\- tt+1 /2 
2 

r(-a+l/2) ^\ Q -V2 
.2> 



(54) 



we obtain 
1 - A 



Substituting (|4]) into > 

pax) « ^-V 2 ) 

tl j r(a) V 27rflt 
where we introduced the coefficient 
|r(-a+ 1/2 



2a-l 



A 



(47) Eq. 



r(a- 1/2) 
^) can be written in the form 

lnPt(ar)-lnP t (0)«-A(|) 



2q-1 



(55) 



(56) 



(57) 



We see that In P t (x) approaches x = as a power of x 
lower than 2 (for 1/2 < a < 3/2). The slope d\nP/dx 
at x — > is zero for a > 1 and infinite for a < 1. 



VIII. COMPARISON WITH THE DOW-JONES 
TIME SERIES 

To test the model against financial data, we down- 
loaded daily closing values of the Dow- Jones industrial 
index for the period of 20 years from 1 January 1982 to 
31 December 2001 from the Web site of Yahoo ||. The 
data set contains 5049 points, which form the time series 
{<S r }, where the integer time variable r is the trading day 
number. We do not filter the data for short days, such 
as those before holidays. 

Given {S T }, we use the following procedure to extract 
the probability density P^ DJ \r) of log-return r for a 
given time lag t. For the fixed t, we calculate the set 
of log-returns {r T = In S T +t/ S T } for all possible times 
r. Then we partition the r-axis into equally spaced bins 
of the width Ar and count the number of log-returns r T 
belonging to each bin. In this process, we omit the bins 
with occupation numbers less than five, because we con- 
sider such a small statistics unreliable. Only less than 
1% of the entire data set is omitted in this procedure. 
Dividing the occupation number of each bin by Ar and 
by the total occupation number of all bins, we obtain the 
probability density pj~ DJ \r) for a given time lag t. To 

find p£ D (x), we replace r — > x + pt. 

Assuming that the system is ergodic, so that ensem- 
ble averaging is equivalent to time averaging, we com- 
pare p( DJ ^ (x) extracted from the time-series data and 
Pt(x) calculated in previous sections, which describes en- 
semble distribution. In the language of mathematical 
statistics, we compare our theoretically derived popula- 
tion distribution with the sample distribution extracted 
from the time series data. We determine parameters 
of the model by minimizing the mean-square deviation 
E xt \ lnP t DJ \x) - lnP t (a;)| 2 , where the sum is taken 
over all available x and t = 1, 5, 20, 40, and 250 days. 
These values of t are selected because they represent dif- 
ferent regimes: jt <C 1 for t — 1 and 5 days, jt w 1 for 
t = 20 days, and yt > 1 for t — 40 and 250 days. As 
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TABLE I: Parameters of the Heston model obtained from the 
fit of the Dow- Jones data using p = for the correlation coeffi- 
cient. We also find 1/7 = 22.2 trading days for the relaxation 
time of variance, a = 2j9/k, 2 = 1.3 for the parameter in the 
variance distribution function and xo = /t/7 = 5.4% for 
the characteristic scale (A4) of x. 



Dow-Jones data, 1982-2001 



Units 



1/day 



1 / year 



7 



4.50 x 10~ 2 8.62 x 10~ 5 2.45 x 10~ 3 5.67 x 10' 



11.35 



0.022 



0.618 



0.143 



Figs, g and || illustrate, our expression (J^) and ( p9| ) for 
the probability density Pt(x) agrees with the data very 
well, not only for the selected five values of time t, but 
for the whole time interval from 1 to 250 trading days. 
However, we cannot extend this comparison to t longer 
than 250 days, which is approximately 1/20 of the entire 
range of the data set, because we cannot reliably extract 
P^ DJ \x) from the data when t is too long. 

The values obtained for the four fitting parameters (7, 
9, k, p) are given in Table |. We find that our fits are 
not very sensitive to the value of p, so we cannot reliably 
determine it. Thus, we use p — for simplicity, which 
gives a good fit of the data. On the other hand, a nonzero 
value of p was found in p5[ | by fitting the leverage corre- 
lation function introduced in 26 and in ju|, |lj, by 
fitting the option prices. 

All four parameters (7, 9, n, p) shown in Table | have 
the dimensionality of 1/time. The first line of the Table 
gives their values in the units of 1/day, as originally de- 
termined in our fit. The second line shows the annualized 
values of the parameters in the units of 1/year, where we 
utilize the average number of 252.5 trading days per cal- 
endar year to make the conversion. The relaxation time 
of variance is equal to I/7 = 22.2 trading days = 4.4 
weeks ~ 1 month, where we took into account that 1 
week — 5 trading days. Thus, we find that variance has 
a rather long relaxation time, of the order of one month, 
which is in agreement with the conclusion of Ref . |2^| . 

The effective growth rate of stock prices is determined 
by the coordinate r m (t) where the probability density 
Pt(fm) is maximal. Using the relation r m = x, m + pt 
and Eq. we find that the actual growth rate is 

p = p — 7^/2^0 ~ p — 9/2 — 13% per year. [Here we 
took into account that luq ~ 7, because 7 3> k/2 in Eq. 
(|32|).] This number coincides with the average growth 
rate of the Dow- Jones index obtained by a simple fit of 
the time series {S T } with an exponential function of t, as 
shown in Fig. (?]. The effective stock growth rate p is com- 
parable with the average stock volatility after one year 
fj = \/9 = 14.7%. Moreover, as Fig. [I] shows, the distri- 
bution of variance is broad, and the variation of variance 
is comparable to its average value 9. Thus, even though 
the average growth rate of stock index is positive, there 
is a substantial probability J_ dr Pt(r) = 17.7% to have 
negative return for t = 1 year. 

According to (f45|), the asymmetry between the slopes 




90 92 94 
Year 



FIG. 7: Time dependence of the Dow- Jones index shown in 
the log-linear scale. The straight line represents the average 
exponential growth in time. 



of exponential tails for positive and negative x is given 
by the parameter pn, which is equal to 1/2 when p = 
(see also the discussion of Eq. (Al) in Appendix |A|). 
The origin of this asymmetry can be traced back to the 
transformation from ([!]) to ([|) using Ito's formula. It 
produces the term 0.5vt dt in the r.h.s. of (0), which then 
generates the second term in the r.h.s. of (g). The latter 
term is the only source of asymmetry in x of Pt{x) when 
p = 0. However, in practice, the asymmetry of the slopes 
Po = 1/2 is quite small (about 2.7%) compared with the 
average slope qf uj /k l/x = 18.4. 

By fitting the Dow-Jones data to our formula, we im- 
plicitly assumed that the parameters of the stochastic 
process (7, 9, k, p) do not change in time. While this as- 
sumption may be reasonable for a limited time interval, 
the parameters generally could change in time. The time 
interval of our fit, 1982-2001, includes the crash of 1987, 
so one might expect that the parameters of the fit would 
change if we use a different interval. To verify this conjec- 
ture, in Figs. and||, we also compare the data points for 
the time interval 1990-2001 with the theoretical curves 
produced using the same values for the parameters as 
shown in Table |. Although the empirical data points 
in the tails for long time lags decrease somewhat faster 
than the theory predicts, the overall agreement is quite 
reasonable. We find that changing the values of the fit- 
ting parameters does not visibly improve the agreement. 
Thus, we conclude that the parameters of the Heston 
stochastic process are essentially the same for 1980s and 
1990s. Apparently, the crash of 1987 produced little ef- 
fect on the probability distribution of returns, because 
the stock market quickly resumed its overall growth. On 
the other hand, the study |2^] indicates that the data for 
2000s do not follow our theoretical curves with the same 
fitting parameters. The main difference appears to be 
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in the average growth rate /i, which became negative in 
2000s, as opposed to +13% per year in 1980s and 1990s. 
Unfortunately, the statistics for 2000s is limited, because 
we have only few years. Nevertheless, it does seem to 
indicate that in 2000s the stock market switched to a 
different regime compared with 1980s and 1990s. 



IX. DISCUSSION AND CONCLUSIONS 

We derived an analytical solution for the PDF Pt(x) 
of log-returns X clS 0j function of time t for the Heston 
model of a geometrical Brownian motion with stochas- 
tic variance. The final result has the form of a one- 
dimensional Fourier integral (|2|) and (H|). (In the case 
p = 0, the equations have the simpler form presented 
in Appendix |A|.) Our result agrees very well with the 
Dow- Jones data, as shown in Fig. |[ Comparing the 
theory and the data, we determine the four (non-zero) 
fitting parameters of the model, particularly the vari- 
ance relaxation time 1/7 = 22.2 days. For time longer 
than 1/7, our theory predicts the scaling behavior ( pq ) 
and ( |37| ) , which the Dow- Jones data indeed exhibits over 
seven orders of magnitude, as shown in Fig. |^. The scal- 
ing function P*(z) = K\{z)/z is expressed in terms of 
the first-order modified Bessel function K\ . Previous es- 
timates in literature of the relaxation time of volatility 
using various indirect indicators range from 1.5 days ||, 
p. 80] to 73 days for the half-life of the Dow-Jones index 
. Since we have a good fit of the entire family of PDFs 
for time lags from 1 to 250 trading days, we believe that 
our estimate, 22.2 days, is more reliable. A close value 
of 19.6 days was found in Ref. [ p5| . 

An alternative point of view in literature is that the 
time evolution of volatility is not characterized by a sin- 
gle relaxation rate. As shown in Appendix the vari- 
ance correlation function Cj ( |C4[ ) in the Heston model 
has a simple exponential decay in time. However, the 
analysis of financial data ^ p. 70] indicates that the 
correlation function has a power-law dependence or su- 
perposition of two (or more) exponentials with the re- 
laxation times of less than one day and more than few 
tens of days. (Large amount of noise in the data makes 
it difficult to give a precise statement.) Ref. |^8| argues 
that volatility relaxation is multifractal and has no char- 
acteristic time. However, one should keep in mind that 
the total range ( p2| ) of variation of is only about 
77% of its saturation value, not many orders of magni- 
tude. As Figs. 2 of Refs. [£9| and [£8| shows, the main 

(v) 

drop of C t takes place within a reasonably well-defined 
and relatively short time, whereas residual relaxation is 
stretched over a very long time. In this situation, a sim- 
ple exponential time dependence, while not exact, may 
account for the main part of relaxation and give a rea- 
sonable approximation for the purposes of our study. Al- 
ternatively, it is possible to generalize the Heston model 
by incorporating more than one relaxation times UM- 



As Fig. g shows, the probability distribution Pt(x) is 
exponential in x for large |a;|, where it is characterized 
by the time-dependent slopes d\nP/dx. The theoretical 
analysis presented in Sec. VI shows that the slopes are de- 
termined by the singularities of the function F t (p x ) from 
( p9| ) in the complex plane of p x that are closest to the 
real axis. The calculated time dependence of the slopes 
d In P/dx, shown in Fig. ^, agrees with the data very well, 
which further supports our statement that 1/7 = 22.2 
days. Exponential tails in the probability distribution of 
stock log-returns have been noticed in literature before 
|| p. 61], however time dependence of the slopes 
has not been recognized and analyzed theoretically. As 
shown in Fig. ^, our equations give the parabolic depen- 
dence of lnP t (x) on x for small x and linear dependence 
for large x, in agreement with the data. Qualitatively 
similar results were found in Ref. |3l|] for a different 
model with stochastic volatility and in agreement with 
the NYSE index daily data. It suggests that the linear 
and parabolic behavior is a generic feature of the mod- 
els with stochastic volatility. In Ref. H, the power-law 
dependence on x of the tails of Pt(x) was emphasized. 
However, the data for S&P 500 were analyzed in Ref. || 
only for short time lags t, typically shorter than one day. 
On the other hand, our data analysis is performed for 
the time lags longer than one day, so the results cannot 
be directly compared. 



Deriving Pt(x) in Sec. IV, we assumed that variance v 
has the stationary gamma- distribution II* (u) (|). This 
assumption should be compared with the data. There 
were numerous attempts in literature to reconstruct the 
probability distribution of volatility from the time-series 
data J32|, [53| . Generally, these papers agree that the cen- 
tral part of the distribution is well described by a log- 
normal distribution, but opinions vary on the fitting of 
the tails. Particularly, Ref. (33| performed a fit with an 
alternative probability distribution of volatility described 
in Ref. p. 88]. Unfortunately, none of these papers at- 
tempted to fit the data using Eq. (|l^) , so we do not have 
a quantitative comparison. Taking into account that we 
only need the integral (p7|), the exact shape of Tl*(v) may 
be not so important, and Eq. (^|) may give a reasonably 
good approximation for our purposes, even if it does not 
fit the tails very precisely. 

Although we tested our model for the Dow- Jones in- 
dex, there is nothing specific in the model which indicates 
that it applies only to stock market data. It would be 
interesting to see how the model performs when applied 
to other time series, for example, the foreign exchange 
data 1 34 , which also seem to exhibit exponential tails. 
The study |27| indicates that our Pt(x) also works very 
well for the S&P 500 and Nasdaq indices for 1980s and 
1990s. However, in the 2000s the average growth rate /j, 
of the stock market changed to a negative value, which 
complicates separation of fluctuations from the overall 
trend. 
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APPENDIX A: THE CASE p = 



Let us expand the integral in (Al) for small x: 



(Bl) 



where the coefficients are the first and the second mo- 
ments of exp[F t (p x )] 



As explained in Sec. VIII , we fit the data using p = 
for simplicity. In this case, by shifting the variable of 
integration in (|28|) p x — > p x + i/2, we find 



P t (x) = e- x ' 2 
where a = 2j6/k 2 , 

F t (p x ) = - a In 
and 



dp x 
2tt 



(Al) 



~2~ 



cosh ■ 



at 



a 2 + 1 2 . , at 

— — sinn — 

2ja 2 



(A2) 



a = V7 2 + k 2 (pS + 1/4) « 7V / 1+pIR7t 2 T- (A3) 

Now the function F t (p x ) is real and symmetric: F t {p x ) = 
F t (~p x ). Thus, the integral in (|Al|) is a symmetric func- 
tion of x, and the only source of asymmetry of Pt(x) in x 
is the exponential prefactor in (|Al|), as discussed at the 



end of Sec. VIII 



In the second equation (A3), we took into account that, 
according to values shown in Table |, n 2 /Aj 2 <C 1. In- 
troducing the dimensionless variables 

i=jt, x = x/x , p x =p x x , x = At/7, (A4) 



Eqs. (Al), (A2), and (A3) can be rewritten as follows 

e -x/2 f+°°dp x 



•t'o 



2tt 



where a = y/l + p\ and 
F~ t (p x ) 



at 
~2 



a In 



cosh ■ 



at 



C ,ip x x+Fj{p. t ) 



a 2 + 1 , , fli 

= — smh — 

2a 2 



(A5) 



(A6) 



It is clear from (A4), ( A5), and (A6) that the parameter 
a determines the shape of the function Pt(x), whereas 
I/7 and xo set the scales of t and x. 

In the limit i 2, the scaling function (|3^) for p = 
can be written as 



P t ( x ) = N t e- x/2 K 1 (z) /z, z= Vx 2 +t 2 , 



(A7) 



where t = at/ 2 = tdjx\ and Nt = te* /ttxq. Notice 



that Eq. (A7) has only two fitt ing para mete rs, xq and 
9, whereas the general formula ( |A5| ) and (A6) has three 
fitting parameters. As follows from (p^), for i>i and 
1, Pt{ x ) « exp(— |x| / xo) , so 1/xq is the slope of the 
exponential tails in x. 



-\-oo -\-oo 

— OO —OO 

. (B2) 

On the other hand, we know that Pt(x) is Gaussian for 
small x. So, we can write 



P t {x) w ^e^e-" 212 / 2 "", 



(B3) 



with the same coefficients as in (Bl). If we ignore the ex- 
istence of fat tails and extrapolate ( p33| ) to x € (—00, 00), 
the total probability contained in such a Gaussian ex- 
trapolation will be 



+00 



(V, = / dxp e- x/2 ^ 2x2/2fH ' = J^l e Mo/8M2_ ( B4 ) 



Obviously, Gt < 1, because the integral (B4) does not 
take into account the probability contained in the fat 
tails. Thus, the difference 1 — Gt can be taken as a mea- 
sure of how much the actual distribution P t (x) deviates 
from a Gaussian function. 



We calculated the moments (B2) numerically for the 
function F given by ( |A2| ) , then determined the Gaussian 
weight Gt from (B4) and plotted it in Fig. [|as a function 
of time. For t — ► 00, Gt — * 1, i.e. Pt(x) becomes Gaussian 
for very long time lags, which is known in literature [0. 
In the opposite limit t — * 0, F t (p x ) becomes a very broad 
function of p x , so we cannot calculate the moments po 
and pi numerically. Th e singular limit t — > is studied 
analytically in Sec. VII. 



APPENDIX C: CORRELATION FUNCTION OF 
VARIANCE 

The correlation function of variance is defined as 



= (v t+T v T ) = / dvi / dvvU t (v\vi)viIl*(vi). 
Jo Jo 

(CI) 

It depends only on the relative time t and does not de- 
pend on the initial time r. The averaging (. . .) is per- 
formed over the ensemble probability distribution, as 
written in (CI), or over the initial time r for time-series 
data. Eq. (CI) has the same structure as in the influence- 
functional formalism of Feynman and Vernon pa] , where 
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nt(w | Vi) represents the conditional probability propaga- 
tor from the initial value Vi to the final value v over the 
time t, and Tl„(vi) represents the stationary, equilibrium 
probability dis trib ution of Uj. 

Using Eqs. (CI) and (g), it is easy to find the limiting 
values of C, 



(C2) 



- = r(l + 0.77), 
a ,' 



where we used the numerical value from Table D. 

Differentiating Eq. (CI) with respect to t and using Eq. 



(^J), we find that C\ v ^ satisfies the following differential 
equation: 



^- = -7(c«-e*). 



(03) 



Thus, C t changes in time exponentially with the relax- 
ation rate 7: 



d v) = e 2 ( 1 



= -7* 



(C4) 
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